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Abstract. Mixing is relevant to many areas of science and engineering, including 
the pharmaceutical and food industries, oceanography, atmospheric sciences, and civil 
engineering. In all these situations one goal is to quantify and often then to improve 
the degree of homogenisation of a substance being stirred, referred to as a passive 
scalar or tracer. A classical measure of mixing is the variance of the concentration of 
the scalar, which can be related to the norm of the concentration field. Recently 
other norms have been used to quantify mixing, in particular the mix-norm as well as 
negative Sobolev norms. These norms have the advantage that unlike variance they 
decay even in the absence of diffusion, and their decay corresponds to the flow being 
mixing in the sense of ergodic theory. General Sobolev norms weigh scalar gradients 
differently, and are known as m,ultiscale norms for mixing. We review the applications 
of such norms to mixing and transport, and show how they can be used to optimise the 
stirring and mixing of a decaying passive scalar. We then review recent work on the 
less-studied case of a continuously-replenished scalar field — the source-sink problem. 
In that case the flows that optimally reduce the norms are associated with transport 
rather than mixing: they push sources onto sinks, and vice versa. 



PACS numbers: 47.51.+a, 47.52.-|-j 
1. Introduction 

One of the most vexing questions about fluid mixing is how to measure it. People 
typically know it when they see it, but specific applications require customised measures. 
For example, a measure might be too fine-grained for some applications that don't 
require thorough mixing. Some measures, such as residence time distributions, are 
designed for open-flow situations where fluid particles are only mixed for a certain 
amount of time. Others, such as the rigorous deflnition of a mixing flow in ergodic 
theory, are better suited to an idealised mathematical treatment. Finally, one of the 
main points of this review is that measures used to quantify mixing in the initial- 
value decaying problem must be interpreted very differently when sources and sinks are 
present. 
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One of the earliest attempts to quantify mixing was by the chemical engineer and 
bomb disposal officer Peter V. Danckwerts [26]. Danckwerts realised that scale was an 
important consideration; he identified the large-scale breakup of fiuid into clumps and 
the subsequent homogenisation at small scales due to diffusion as separate processes [25] : 

The breaking-up and the interdiffusion are, in the case of liquids, largely 
independent processes which produce distinguishable results. The former 
reduces the size of the clumps, while the latter tends to obliterate differences of 
concentration between neighbouring regions of the mixture. It therefore seems 
desirable to use two quantities to describe the degree of mixing — namely the 
scale of segregation and the intensity of segregation. 

Two other pioneers are the oceanographers Carl Eckart [32] and Pierre Welander [92], 
who also identified the complementary roles of mechanical stirring and diffusion. 
Following Eckart, modern parlance refers to these two stages as stirring and mixing, 
the distinguishing feature being that stirring is a mechanical action, whilst mixing is 
diffusion-driven or the result of coarse-graining. Welander in particular was emphatic 
about the important role of stirring in creating filaments, which can subsequently be 
smoothed by diffusion. It was then Batchelor [12] who identified the length scale of the 
filaments, at which stirring and diffusion achieve a balance, now known as the Batchelor 
scale. 

Let us pin a mathematical meaning on the ideas above. Danckwert's scale of 
segregation is a correlation length of the concentration of a mixture. His intensity of 
segregation is a normalised variance of the concentration. It is through the variance that 
the connection to norms first appears, in this case as the L^-norm of the concentration 
field. If 6{x, t) is the concentration of a passive scalar — such as temperature, dye, or 
salt — then the variance is 



are the L^-norm and mean value of 9, respectively, with Vl the spatial domain and \Vt\ 
its volume (or area in two dimensions). 

Why is the variance a good measure of mixing quality? A first answer is that 
it measures fiuctuations from the mean, and a mixed state is exactly one where the 
concentration is equal to the mean — i.e., it is uniform. But there is a second, 
more intimate reason why variance is important: from the classical advection-diffusion 
equation for an incompressible velocity field u{x,t) and diffusion coefficient 



Var^= ll^ll^ 




where 




(1.2) 



dO 
di 



+ u-ve = kAo, 



V • w = 0, 




we find that the concentration variance is monotonically driven to zero in the absence of 
sources (see section 3). Indeed, after a few integrations by parts and assuming boundary 
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conditions that conserve the total amount of 6 (no-flux or periodic, see section 2), we 
find the L^-norm and mean (1.2) obey 

|W=0, |||()f = -2.||vC, (1.4) 

SO that 

^Var^ = -2k||V^|P. (1.5) 

dt " " ^ ' 

Observe that the right-hand side is negative-definite unless 9 — const., i.e., the 
concentration is uniform. Hence, the advection-diffusion equation says that the 
concentration is driven to a uniform state with Var6' = 0, at a rate dictated by the 
product 2/t||V^^||2- Determining this rate is an important aspect of the scalar decay or 
initial value problem, where we have some initial concentration field 6{x,t) = 6q{x) and 
want to know how fast it is mixed by a velocity field u{x,t). There is a vast literature 
focused on determining and estimating this decay rate. Note that the advecting velocity 
field does not appear directly in (1.5): its role is to increase gradients of concentration. 

Thus, monitoring variance is a simple way of quantifying the effectiveness of a 
mixing process. It has, as mentioned above, the dual advantages of being intuitive and 
of being mathematically sound. In what situations, then, is the variance a less-than- 
ideal measure of mixing quality? The answer is: when the difFusivity k is very small 
and the stirring process is very effective. The typical situation in that case is that the 
decay term 2/s:|| V^^||2 becomes independent of k. The physical picture is that the stirring 
sharpens gradients of 6 until they are large enough that diffusion easily smooths them 
out. A balance between advection and diffusion is then reached, and the variance decays 
at its optimal rate. 

In theory, this is fine; in practice, it is disastrous. If we are trying to optimise the 
mixing process, it means we have to keep track of scales down to lengths of order k^/^, 
which can be very small. Our simulations are limited by the small quantity k, and yet 
the decay rate is independent of k,. However, if we try and avoid this by setting k, = 
in (1.3), we find that variance is exactly conserved. The limit is singular: taking k to 
zero is completely different from setting k = 0. 

This problem is similar to what happens for turbulence: letting the viscosity tend 
to zero in the Navier-Stokes equations does not necessarily recover solutions to Euler's 
equations. However, here the problem is less severe, since the concentration 9 does not 
feed back on the flow. Hence, it is sensible to solve the pure advection equation 

B9 

+ u-V9 = (1.6) 

ot 

and try to extract some measure of mixing from 9; it's just not possible to use variance 
as a measure. The key is to use a different norm that downplays the role of small 
scales. For example, a simple choice is || V^^^^||2. The operator V^^ will be defined more 
carefully later, but for now it suffices to understand that it 'smooths out' 6, so that 
small-scale variations are not detected. (We make sure the operator is well-defined by 
restricting to functions with vanishing mean, (9) = 0.) 
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Figure 1. Satellite observations showing concentrations of carbon monoxide (CO) in 
the atmosphere. Red corresponds to high levels of CO (450 parts per billion) and blue 
to low levels (50 ppb). Note the immense clouds due to grassland and forest fires in 
Africa and South America (photo from NASA/NCAR/CSA). 



So far our discussion has been of equation (1.3) when the total quantity of scalar is 
conserved, i.e., {6) = const. This is the case most often discussed in the literature, and 
in recent years much work has gone into understanding the factors affecting the decay 
rate of the concentration and its variance. There is also a vast literature on the case 
where there are sources and sinks, either within the fluid domain or at its boundaries. 
In geophysical and industrial situations this is often more relevant. Figure 1 shows a 
satellite image of the concentration of carbon monoxide in the atmosphere. Observe 
that there are several large and many small sources, distributed in a complex manner 
throughout the globe. This begs the question: what do we mean by well-mixed in a 
case as in figure 1? Unmixed features are persistent, since they are replenished by an 
inhomogeneous source. But we can still use a norm-measure such as the variance: it 
will not tend towards zero with time, but a smaller value of variance still indicates that 
the passive scalar is getting mixed to some degree. A stirring flow that is effective at 
mixing will thus presumably tend to reduce the variance of the concentration field. 

It has become apparent in recent years that the quality of mixing we obtain will 
depend strongly on the source-sink distribution, in addition to depending on the stirring 
flow itself. Let us give a simple but extreme example of this, discussed by Plasting & 
Young [64] and Shaw et al [71]. Consider a two-dimensional biperiodic square domain 
where sources and sinks are present, as in figure 2(a). The equation to be solved is the 
analogue of (1.3) with a body source term: 

do 

— + u- Ve = kA9 + s{x,t), (1.7) 

with periodic boundary conditions. Here we choose a simple source-sink distribution, 
s{x) = sin{2iTx/L). If we fix the kinetic energy of the flow, | UpwH^,:}: what is the 

I We consider incompressible flows with constant density throughout, so that fixing the kinetic energy 
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Figure 2. (a) A source-sink distribution in a two-dimensional periodic square domain, 
(b) Sketch of the velocity field that most effectively reduces the concentration variance 
of the source. 

incompressible stirring velocity field that will most effectively reduce the variance? The 
answer sketched in figure 2(b) is somewhat surprising: the most effective stirring consists 
of a uniform (constant) flow. This flow carries the hot fluid onto the cold sink, and cold 
onto hot, thereby reducing the concentration variance as much as possible. 

The reason this is surprising is that 'common wisdom' in mixing assumes that the 
best stirring is either turbulent or exhibits chaotic trajectories [4,59]. Such complex 
behaviour increases concentration gradients and thus allows diffusion to act more 
effectively. However, this particular source-sink configuration is best mixed not by 
creating small scales, but rather by transporting fluid appropriately. One might not 
call this 'mixing' in the strict sense, but it is dependent on diffusion, as (1.7) does not 
converge to a steady-state without it. 

Mixing a scalar field whose fluctuations are constantly replenished by spatially 
inhomogeneous sources and sinks is a problem with a long history. Townsend [86,87] 
was concerned with the effect of turbulence and molecular diffusion on a line source of 
temperature — a heated filament. The spatial localisation of the source, imposed by 
experimental constraints, enhanced the role of molecular diffusivity. Durbin [31] and 
Drummond [30] introduced stochastic particle models to turbulence modelling, and these 
allowed more detailed studies of the effect of the source on diffusion. Sawford & Hunt [68] 
pointed out that small sources, such as heated filaments, lead to an explicit dependence 
of the variance on molecular diffusivity. Many refinements to these models followed, see 
for instance [2,3,15,85] and the review by Sawford [67]. Chertkov et al [16-20] and 
Balkovsky & Fouxon [8] treated the case of a random, statistically-steady source. 

We now give a brief outline of the review. In section 1.3 we recall some basic 
properties of the advection-diffusion equation, to complement the earlier material in 

is equivalent to fixing the i^-norm For that reason, we shall often refer to ||tt||^ itself as the 

kinetic energy. 
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this introduction. Section 3 is devoted to a review of Sobolev norms. The rest of the 
review is divided into two parts. Part I discusses norms as measures of mixing for the 
freely-decaying passive scalar. Section 4 connects norms with the concept of mixing in 
the sense of ergodic theory. In section 5 we show how negative Sobolev norms can be 
used to optimise flows to achieve rapid mixing. 

Part II, which forms the bulk of the paper, is devoted to advection and diffusion in 
the presence of sources and sinks. In section 6 we introduce mixing efficiencies, measures 
of mixing based on norms. We give some upper bounds on these efficiencies in section 4. 
In section 8 we investigate the dependence of efficiencies on functional features of the 
source-sink distribution. We derive mixing efficiencies from a homogenisation theory 
approach in section 9. In section 10 we discuss optimisation of mixing efficiencies. 
Finally, we offer some closing comments in section 11. 

2. The advection-diffusion equation 

The main equation discussed in this review is the advection-diffusion equation for a 
passive scalar with concentration field 9{x,t), 



Following the usual route, we will assume that the domain of interest f2 is a periodic 
square box with dimension d, with spatial period L in each direction. Of course, 
everything discussed in this paper can be repeated for a more general closed domain 
with no-flux boundary conditions, but this adds little to the discussion. In particular, 
by using a periodic domain we can use Fourier series expansions, which makes many 
calculations explicit. 

The mean {9) satisfies 




(2.1) 




(2.2) 



with solution 




(2.3) 



Thus, if we replace ^ by a new variable 




(2.4) 



then e' 



obeys the modified equation 




(2.5) 



The new concentration field 9' and source s' have spatial mean zero. For the remainder 
of this paper, we drop the primes and assume without loss of generality that both 9 
and s have zero spatial mean. 
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3. Norms 

3.1. Definitions and basic properties 

In this section we introduce the measures of mixing we'll be using for the rest of the 
paper. The Sobolev norm we use for the space H'^{fl) is 




or in terms of Fourier series, 

\\f\\H.-{j2(^ + k'Ly\U\'y\ (3.2) 

k 

where k := |A;| and fk are the Fourier coefficients. In this review we will prefer to use 
the seminorm on the homogeneous space i7^(fi) (note the dot over H), 




or in terms of Fourier series 

fe 

If wc take a function / with Fourier coefficients behaving asymptotically as \fk\ ~ k^, 
/c ^ 1, then the norms (3.1)-(3.4) converge (exist) for g + p < —d/2, where d is the 
dimension of space. 

The norm (3.3) has a more intimate connection with solutions of the advection- 
diffusion equation than (3.1), as will be described later (section 3.2). Note that the 
manner in which we have defined (3.1) and (3.3) allows for q positive, negative, or even 
fractional. In section 2 we showed that we could assume (9) — 0, so we restrict attention 
to functions / with mean zero. In that case (3.3) becomes a true norm since ||/||^q = 
if and only if / is zero. 

In fact it docs not matter much which of the two norms (3.1) or (3.3) we 
use, since they are equivalent for zero-mean functions: by Poincare's inequality, we 
have ||(-A)5/V|| > (W^)^ ll/ll for ? > 0, so that 

11/11^. < L-'^ WfWm < (1 + C^nry/' WfU, , q>0, (3.5) 
for all zero-mean functions /. For q < Poincare's inequality is reversed, so we have 

(1 + {2nrY' Wfh. < L"^ WfWm < WfWn. , ? < 0- (3.6) 

Equivalence means that if one of the two equivalent norms goes to zero, then the other 
must as well, and they must do so at the same rate [54,55]. 
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For the mathematically minded, we can give a rigorous definition of what kind of 
functions live in the negative Sobolcv space i?^, with q < 0, given that we understand 
the space H~'^. The space H'^ is defined as the space dual to with respect to the 
standard pairing 

if^9)-TL f f{x)g{x)dn, (3.7) 

with the dual norm 

||/||^_,. = sup ^4f^ . (3.8) 
geH-1 \\9\\h-i 

The norm (3.3) is equal to the dual norm (3.8). To show this, first observe that by using 
the Cauchy-Schwarz inequality, 

llfi'llii"-'? 

independent of ^, so ||/||^_,, < ||/||^,. Now let g = {-Afif: 
{f,g) ifA-^r^f) ll(-A)Vir 



\gU., ||(-A)2«/||^_, ||(-A)^/|| 



Hi 



(3.10) 



Since the dual norm is defined as a sup over g, we have > ||/||//9- We conclude 

that 11/11^-9. = ll/llij?- The same argument can also be used to show H/Hj:^-,* = WIWhi 
for the inhomogeneous spaces. 

3.2. Evolution in time 

To get a feel for what these norms are telling us about mixing, it is helpful to examine 
how they evolve in time. That is, given that 9 obeys (2.1), what is dWOW^jq/dtl We 
start from 

\eW%,^{(-A)wdte) (3.11) 



1^ II/1II2 
^dt 



where (•) denotes an average over the periodic domain ^2. Inserting (2.1) for dfO, 
^ d 



= - {{-AyOu ■ve)-K ||(-A)(''+^)/2^|f + ((-A)'^^s) . 

The case q — gives the evolution of the variance, for which the velocity term on the 
right integrates away: 

lf^\m%o^-Kmf+{es). (3.12) 

There are two other cases that give particularly nice equations. The case q = I gives 
the evolution of scalar concentration gradients: 

1^ ii^ii^^ ^ _ . vu • V9) - n IIA^II' - {AOs) . (3.13) 
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The first term on the right is the famihar 'stretching' term, which says that gradients 
are increased or decreased proportionally to their alignment with the principal axes of 
the rate-of-strain tensor. A direction of positive strain will decrease gradients, whilst a 
direction of negative strain will increase gradients. 

The final case of interest to us is g = —1, for which 

\^ \\e\\\_, = {V-^ . Vu ■ V-^e) - n ll^ll' - (A-^^s> . (3.14) 

Here we interpret via its action on Fourier modes 

(V-^^)fe = -^4, (3.15) 

so that V ■ V^^^ = 9. (Recall that we are restricting to functions with vanishing mean.) 

Compare the first term on the right-hand side of (3.14) to the same term for (3.13): 
velocity gradients have the opposite effect on ^ ||6'||^_i as they do on ^ || 6*11^^1- This 
is intuitively clear: the creation of concentration gradients will tend to make 6 very 
filamented. We will see in section 4 that this will cause it to converge weakly to zero, 
and that this implies that any negative Sobolev norm must go to zero. Another argument 
that the two norms should evolve with opposite trends arises from 

\\e\\^ ^-(ye- v~^e) < pWh, ||^||^_i . (3.16) 

Since the variance ||6'||^ is conserved when k — s — 0, inequality (3.16) implies that if 
the norm converges to zero, then must diverge. 



Part I: The decaying problem 

4. Mixing in the sense of ergodic theory 

A divergence-free velocity field u{x,t) generates a time-dependent function 6{x,t) via 
the advection equation, 

— + u-Ve^O, V-u = 0, (4.1) 

at 

with initial condition 6' (cc, 0) = 6o{x) G L^{Q). Note that (4.1) preserves ||6'(-,t)|| = ||6'o|| 
for all time, so ||^(-,t)|| is uniformly bounded in time by || 6*011. We have left out diffusion 
in (4.1): in the present section we will discuss how we can define mixing without 
appealing to diffusion, and how the norms (3.1) and (3.3) can be related to this type of 
mixing. The definition of mixing we will use is the one from ergodic theory. 

The property of mixing in the sense of ergodic theory is a little opaque when 
described mathematically. The rigorous definition is as follows: 

Let {X,A,ii) be a normalised measure space and : X ^ X be a measure- 
preserving flow. is called mixing if 

lim ii(A n S-\B)) = ii{A) ii{B), for all A, Be A. 
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Here X is our domain Q, ^1 is a so-called a- algebra over X, and /x is a measure. The 
elements of A are measurable sets, which wc can think of as patches or 'blobs' in Vt. The 
measure /i assigns a positive real number to a set in A. We have been using Lebesgue 
measure — the dfi that appears in our integrals — which means the measure of a blob 
is just its volume. The flow 5* for us takes Oq to 0{-,t) subject to (4.1). The flow 5* 
preserves volume, and hence Lebesgue measure. 

Intuitively, the definition of mixing works as follows. As we only deal with reversible 
systems here, we can replace S~*{B) by S*'{B) in the definition above, since such a system 
should be mixing both forward and backward in time (this makes things conceptually 
easier). Think of A as a fixed reference patch, and B as a. blob that gets stirred. Since 
the transformation is volume-preserving, S^{B) might stretch and filament, but it does 
not change its volume. However, if it fills the domain 'uniformly,' then the volume of 
its intersection with A just ends up being proportional to A and -B's volume. Crucially, 
this is true for every set A and i?, so everything ends up everywhere. (Note that this is 
stronger than ergodicity, which only requires sets to visit every point, but not necessarily 
be 'everywhere at the same time.') 

Mathew et al [54,55] have introduced a norm, called the mix-norm, which captures 
the property of mixing in the sense of ergodic theory. The mix-norm is somewhat 
cumbersome to define on the torus, but to give an idea of its flavor we will describe it 
for the one-dimensional periodic interval [0, L]. First, deflne 



for all X, w G [0,^^]. The function d{9,x,w) is the mean value of the concentration 
9 in an interval of width w centred on x. The mix-norm $(6') is then obtained by 
averaging o?^ over x and w: 



In words, the mix-norm averages the concentration over an interval of width then 
averages the square of this over all intervals and all widths. In dimensions greater than 
one, the deflnition of the mix-norm involves integrals over balls of varying sizes instead 
of intervals. 

For our purposes here it suffices that the mix-norm (4.3) and its higher-dimensional 
generalisation are equivalent to the norm (3.1) with q = —1/2 [55]. We will also see below 
that all Sobolev norms with q <Q capture the property of mixing in the same sense as 
the mix- norm, as shown by Lin et al [50]. For that reason, we extend Mathew et al 's 
terminology and often refer to any negative Sobolev norm as a mix-norm, not just for 
the case q — —1/2. The term multiscale norm encompasses the Sobolev norms for any 
value of g — positive, negative, or zero. 

The connection between negative Sobolev norms and mixing involves the property 
of weak convergence. We refer the reader to the book by Lasota & Mackey [48] for a 




x—w/2 



(4.2) 




(4.3) 
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more complete discussion of the relation between weak convergence and mixing. Weak 
convergence is defined as follows: 

A time- dependent function f{x,t), f{-,t) e L'^{^), is weakly convergent 

to /oo e L\n) if 

\im{f{;t),g) = (/oo,^), for all g e L\Q). 

r— ^-oo 

Note that if f{-,t) converges weakly to /oo, then f{-,t) — /oo converges weakly to zero. 

Instead of a time-dependent function, wc can also use a discrete sequence {cm}, 
where m is like time. A simple example of a sequence that converges weakly to zero is 
given by Cm — sin(27rma;/L), since 

lim / sin(27rmx/L) g(x) dx — 

for all g G L'^{Q), by the Riemann-Lebesgue lemma. The connection to mixing is 
that under stirring a passive scalar usually develops finer and finer scales, much like 
the function sin(27rma;/L) with increasing m. In practice diffusion smooths out these 
large gradients and the concentration field tends to zero at every point. However, if 
we ignore diffusion and retain the small scales we can still detect this mixing process 
by 'projecting' onto test functions such as g. We now discuss the connection between 
mixing in the sense of ergodic theory and the norms (3.3). 

The following is a slightly more general version of the mix-norm theorem by 
Mathew et al [55] (the case g = — 1/2 is equivalent to their theorem): 

A time- dependent function f{x,t), where f{-,t) G L^(0) has mean zero and is 
bounded in the norm uniformly in time, is weakly convergent to zero if and 
only if 

lim ||/(-,t)||j:^, = 0, foranyq<0. 

r— >oo 

The proof given by Lin et al [50] is reproduced in Appendix A. 

Another theorem from Mathew et al [55] now implies that the dynamics generated 
by u{x,t) are mixing in the sense of ergodic theory if and only if limj^oo ||^(', ^)|lii-9 = 0, 
for any q < 0. This is a direct consequence of the mix-norm theorem above. 

Note that the equivalence of the norms q < 0, with mixing in the sense of 

ergodic theory is only a useful concept for the freely-decaying problem. In the presence 
of sources and sinks, diffusion plays the essential role of making an asymptotic state 
possible (in its absence solutions can diverge), so we cannot simply solve (4.1) with a 
source term on the right and expect to get anything sensible. 

We close this section with a rule of thumb to interpret the decay of mix-norms. 
Consider a function / with Fourier coefficients fk, where the coefficients vanish when k 
contains odd wavenumbers. Now define /' by f'f^ = f(k/2), that is, /' is the same as / 
but with all scales divided by two. Prom definition (3.4), we have 

WfWm = (E^'' IhM'Y' = (E(2^')'' l^'l')'^' = 2^ WfWn^ ■ (4-4) 
fe fe' 
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Thus, a refinement of scales by a factor of two leads to a decrease in ||-||^9 of a factor 2'^ 
[q < 0). For q = —1, the norm decreases by half when scales are refined by half. 
For g = -1/2, the mix-variance ||'||^q decreases by half when scales are refined by half. 
In both cases, the norms decrease at a rate which reflects the creation of small scales. 



5. Optimisation for decaying problem 
5.1. Optimal control 

Mathew et al [54] have used optimal control techniques [9,43,45,73] to find velocity 
fields that rapidly reduce the norm H~^/'^ of a concentration field. (See also [10, 23, 24, 
40,51,69,76,91].) The energy of the fiow is held fixed (more precisely, its total action 
over a time interval). They assume that the velocity field can be expressed as a linear 
combination of steady incompressible velocity fields Ui{x) as 

n 

u{x,t) — '^^ai{t)ui{x), V • Uj = 0. (5.1) 

i=l 

The coefficients «j(t) the controls that are adjusted to achieve the optimisation. It is 
assumed that the flow (5.1) can be realised in practice for a given set of functions a,i{t). 
The quantity to be optimised is the Sobolev norm ||^^||j:^-i/2 of a concentration fleld 
satisfying the advcction equation (1.6), for some initial condition 9q{x). 

In this formalism, the time-integrated energy (i.e., the action A) and the advection 
equation (1.6) itself enter an augmented functional as constraints: 

W[oc, e, z, 7^] = \\e{;tf) -z(^A- J\{t)-R- cx{t) dt^ 



where z and r]{x,t) are Lagrange multipliers, and the matrix R describes the kinetic 
energy for the individual velocity fields in (5.1), 

Rij := I {ui{x) Uj{x)) . (5.3) 

Note that the functional (5.2) involves the norm of the concentration field only at the 
final time tf. However, the constraints involve the entire history of 9{x,t). 

Taking the functional (Prechet) derivatives of (5.2) with respect to a, 9, z, rj and 
equating to zero leads to a two-point boundary value problem: 9{x,t) is specified at 
the initial time t = 0, but r){x,t) is specified at the final time tf (sec Eq. (3.2) in [54]). 
Both 9 and 77 (the 'costate field') satisfy advection equations of the form (1.6), so they 
can be solved by following particles on Lagrangian trajectories (backwards in time for 77), 
following the velocity field given by the current best guess for a{t). This guess can be 
varied following an iterative procedure, such as the conjugate gradient method. 
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? = 0.6 / = 0.8 r=1.0 



Figure 3. Snapshots of the concentration field 6{x^ t) advected by the optimal solution 
for the velocity fields defined by (5.4) (from Mathew et al [54]). 

As an illustration of the method, Mathew et al [54] apply it to the velocity fields Ui 
and U2 with streamfunctions 

'?/'i(x, ?/) = sinx siny, ■?/'2(a;, = cos x cosy. (5.4) 

They set the time interval tf = 1 and action A = 1/4, with initial concentration 6q{x) = 
sin y. Figure 3 shows the concentration field evolved by their optimal solution, which is 
a time-dependent linear combination of (5.4). Note that individually the two fiows (5.4) 
lead to poor mixing, since they are steady and two-dimensional. Figure 3 exhibits finer 
and finer scales as time evolves, a hallmark of chaotic advection, as well as a roughly 
exponential decay of the Sobolev norm ||^^||j:^-i/2 (Figure 5.2), dotted line). 

5.2. Local-in-time optimisation 

Lin et al [50] proposed an alternative to the global optimal control approach of 
Mathew et al [54]. Instead of focusing on the norm of the concentration field, they 

instead examined H^^. As the theorem in section 4 tells us, any negative Sobolev norm 
will capture mixing in the sense of ergodic theory, so there is no profound difference in 
using either norm. However, the rate at which a different norms decrease will in general 
be different, so different optimal solutions can be obtained. The advantage of arises 
from the evolution equation (3.14), in the absence of sources and diffusion: 

^||^||i,. = 2(V-^^-V«-V-^^>. (5.5) 



Using multiscale norms to quantify mixing and transport 



14 



The right-hand side is a simple expression that can easily be extremised instantaneously ^ 
in the sense that given 9 at any instant we tweak the velocity to cause the norm to decay 
as fast as possible. This local-in-time approach can never do better than global optimal 
control, but it is often good enough, as we will see. But most importantly, it is much 
less computationally expensive since we do not have to 'peer into the future' and evolve 
the system forward in time to determine the optimal current velocity field. 

In order to formulate an optimisation problem, we must impose some constraints 
on the velocity field. In section 5.1 we imposed fixed total kinetic energy through (5.3). 
Now we will consider two types of constraints, fixed energy or fixed power. For a 
Newtonian fluid, the power is proportional to the norm ||Vw||^, but we shall refer 
to this integral as 'the power' even if the fluid is not of this type. The two constraints, 
then, are to respectively fix 

(fixed energy) (5.6) 



or 



d 

\\^uf = J2 {i^iUjf) = (fixed power). (5.7) 

These define the root-mean-square velocity U and rate of strain 7"^ of the stirring. 

We now proceed with the optimisation technique, that is, to maximise the right- 
hand side of (5.5). With a few integrations by parts we recast (5.5) in the form 

^ \\e\\l_, = -2{9u- V{A'^e)) = -2{u- F{9V(p)) (5.8) 

where is the filtered scalar field, 

cP{x,t):={A-'e){x,t), (5.9) 

and P(-) is the projector onto divergence-free fields defined by 

P(v) :=-y- VA-^(V-'y). (5.10) 

Then with either the fixed energy (5.6) or fixed power (5.7) constraint the velocity field 
maximising the decay rate of is 

Ue = U i,^^^^*^^,, (fixed energy) (5.11) 

||P(^V0)|| \ K J 

or 

^ "7 ||P(^V0)||,-. ^^"^^ ^'-^'^ 

as long as the denominator does not vanish. Hence, or Up is the best stirring velocity 
fields to use at any instant in time, unless the denominator vanishes. However, if either 
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of the norms in the denominators vanishes then P(^V0) = throughout the domain 
and no incompressible flow can instantaneously decrease the norm. For example, 
this will happen if the concentration field satisfies = F{9), which includes cases 
where 9 is an eigenfunction of the Laplacian. If this situation arises in the course of the 
time-evolution of 6*, then some other optimisation strategy must be adopted. 

The most natural alternative when P(^V0) = is to carry to optimisation to the 
order, that is, find the velocity field that minimises 

— \\e\\\.,^2{[u-V(t)Ve-u-{u-ve)/^-\u-ve)\). (5.13) 

(Jj V 

Then the optimal incompressible flow u solves the eigenvalue problem 

Au = P(^(u • V^)V0 + {u ■ V0)V^ - 2[A-\u ■ V9)]Ve^ (5.14) 

for the fixed energy constraint (5.6) or 

A u = -A-^P(^(w • V^) V0 + (u ■ V0) - 2[A-i(w • V^)] V^) (5.15) 

for the fixed power constraint (5.6). In either case we seek the eigenfunction 
corresponding to the minimum eigenvalue A_ < to use as the stirring field 
momentarily, until P(^V0) ^ 0. The eigenvalue problems in (5.14) and (5.15) are 
generally difficult; see [50] for a discussion. In practice, we may need to solve one of 
the eigenvalue problems if we choose a 'bad' initial condition (such as an eigenfunction 
of the Laplacian, which is commonly done), but one it has started the optimisation 
procedure does not seem to get stuck very often. 

This local-in-time optimal stirring strategy is a limiting case of the short-horizon 
optimisation studied by Cortelezzi et al [23] when the horizon becomes infinitesimal, 
but the locality and simplicity allows a much broader class of flows to be used. In order 
to implement it in practice the full scalar field must be monitored so that the optimal 
flow field can be computed at each instant. 

We reproduce here the tests of this optimal stirring strategy presented in Lin et 
al [50]. They used initial scalar distribution ^o(^) = sinx in a domain of size L = 27r in 
d = 2 spatial dimensions, for the fixed power constraint (5.7) with 7^^ = 6.25 x (27r)^, 
equivalent to the amplitude of the bi-component control used by Mathcw et al [54]. 
The results for various norms are shown in Figure 5.2. The optimisation was performed 
for H~^, but the H"^/'^ norm is also plotted to allow a direct comparison with Mathew et 
al [54] . The local-time-optimisation seems to outperform the global optimal control, but 
this is because the former has access to all possible incompressible velocity fields. The 
difference is evident when comparing figure 3 and figure 5, which shows the concentration 
field. The Lin et al [50] solution uses much smaller velocity scales (though always at 
fixed power, so the fiow must slow down). Note also the the optimal fiow in figure 5 is 
suggestively self-similar in time. 




Figure 4. Evolution of norms with the fixed power constraint (5.7) for dQ{x) = s'mx. 
All norms are rescaled by their initial values, and the conserved norm is monitored 
as a numerical check. The optimisation is over all possible velocity fields satisfying the 
power constraint, which is why the local-in-timc optimisation outperforms the optimal 
control approach of Mathew et al [54]. Snapshots of the velocity field are shown in 
figure 5, and in figure 3 for the Mathew et al solution (from Lin et al [50]). 



t = t = 0.2 t = QA 




Figure 5. Snapshots of the scalar field in [0, 2ttY for the fixed-power optimal mixing 
strategy (5.12), with initial condition 9o{x) = sin a; (from Lin et al [50]). 
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6. Mixing efficiencies 

There are several ways to 'calibrate' the norms to measure some efficiency of mixing. 
For the freely-decaying case (no sources or sinks), one can normalise the norm by its 
initial condition to obtain the ratio / H^oIIh?- Foi' the case g = this is Danckerts' 

'intensity of segregation' [25]. The goal of optimisation is then to reduce / H^^ollij? 

as rapidly as possible. This is a strategy that has been used by many authors for the 
case of variance {q = 0) [22] and for negative Sobolev norms [40,50,54], as we discussed 
in section 5. 

In the presence of sources and sinks, the norms ||^||^q typically reach an asymptotic 
steady state (or at least a statistically-steady state) . In that case normahsing the norms 
by their initial values is not helpful, since the asymptotic state is usually independent 
of initial condition. Instead, a convenient measure of mixing efficiency is to normalise 
the time-asymptotic norm by the value it would have in the absence of stirring. 

We define mixing efficiencies (or mixing enhancement factors) by 

(6.1) 

where 9 is the steady solution to (2.1) and 9 is the steady solution to the diffusion 
equation 

(W ~ 

— ^kA9 + s. (6.2) 
ot 

(If the velocity field or source are explicitly time-dependent, then an appropriate long- 
time average must be added to the norms in (6.1); we will see this in section 7.) 

The efficiencies measure the amount by which a norm is decreased by stirring. 
If stirring decreases a norm ||^||^q over its purely-diffusive value, then Eg is larger. 
For q < 0, an increase in efficiency is associated with better mixing, since the flow has 
suppressed fluctuations. 

We might expect that stirring should always decrease the norms from their purely- 
diffusive value. For £i, this is easily shown to be the case [71]. From the deflnition (3.15) 
of the inverse gradient of a mean-zero function, we have 

K \\9\\%^ = {9s) = V • V-'s> = - (Ve • V-'s> 

< ll^ll^i ||V-^s|| (6.3) 

where we used the Cauchy-Schwarz inequality. The steady-state solution of (6.2) 
is ^= -K-^A-h, so V9 = -K-^Vh. We conclude from (6.3) that < or 

£i > 1. (6.4) 

Thus, the efficiency deflned with gradients of 9 is always decreased by stirring. This is 
somewhat counter-intuitive, since we expect stirring to create gradients, but this result 
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holds only for the steady state (or a long-time average). In fact the conventional wisdom 
in mixing holds that stirring creates sharp gradients, and that those sharp gradients 
are responsible for good mixing (see introduction). The bound (6.4) shows that this 
viewpoint must be qualified when sources and sinks are present: stirring may indeed 
create small scales, but the gradients are never as sharp overall as those that would 
build up if we didn't stir at all. 

Perhaps even more surprising is that the efficiencies £o and £_i are not always 
increased by stirring. The possibility of this was mentioned by Shaw et al [71] . Indeed, 
the combination of flow and source-sink distribution given by§ 

u = (sin 2a; cos 2y , — cos 2x sin 2y) , (6.5a) 
s = cos 2x sin y, (6.5b) 

with Pc = 10 has Sq — -991, £_i ~ .638, both less than unity. Here the Peclet number Pe 
is defined in terms of the norm of u as 

Pe := UL/k, U = \\u\\ , (6.6) 

and L — 27r, U — l/-\/2 for the velocity field (6.5a). (These unmixing fiows may be 
related to fiows that create 'hotspots' [36,42].) 

We can optimise the stirring velocity field to give the worst possible mixing 
efficiency for the source (6.5b), using the same techniques as in section 10.2. Figure 6 
shows the resulting streamfunction for Pe = 10, as well as the solutions 9 and 9 to the 
advection-diffusion and diffusion equation, respectively. The optimised unmixing flow 
has £o — .945, £_i ~ .642, a modest improvement over (6.5a) (in fact £_i went up, 
since the flow was optimised for smallest £o). At larger Peclet number the optimised 
solution has lower suboptimal efficiency. The fiow appears to achieve this low efficiency 
by stretching parts of the source along the vertical direction, creating thin peaks, but 
avoiding concentrating it in other places. Such 'unmixing fiows' are fairly rare and 
delicate to construct. It remains true that for a given source most velocity fields will 
have efficiency greater than one, though the specific details of this question have not 
been thoroughly investigated. 

In fact it was thought that sources of the form (6.5b), where the source is an 
eigenfunction of the Laplacian operator, could not lead to £o or £_i less than unity. 
However, the argument in [71] contains a fiaw. To obtain a lower bound on £o, the 
authors solved the constrained optimisation problem 

{9^) < max {(i?') I k{\V^^) = {s^)} (6.7) 

using the Euler-Lagrange equation 

2i?* + 2nnM^ + //s = (6.8) 

§ The question of whether such 'unmixing' flows exist was posed by Charles R. Doering at the 

Workshop on Transport and Mixing in Complex and Turbulent Flows, Institute for Mathematics and 
its AppUcations, MinneapoHs, in April 2010. The form (6.5) is derived from a solution suggested by 
Jeffrey B. Weiss by the end of the workshop. 
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Figure 6. (a) Streamlines of the optimised unmixing flow for the source-sink 
distribution (6.5b) (shown in background) at Pe = 10. (b)-(c) Corresponding 
temperature field with and without stirring, respectively. 



where is the Lagrange multipher enforcing the constraint in (6.7). In terms of the 
Fourier coefficients the solution of (6.8) is straightforward, 



'&*k = IfiSk/if^Kk'^ - 1) 
but /i is the solution of the generally-difficult problem 

2 - iiKk"^ 



11 



Pfcl 



0. 



However, a convexity argument shows that (6.9)-(6.10) is only a maximum if 

where K = 2tt/L is the magnitude of the smallest wavenumber. 

When the source is an eigenfunction of the Laplacian, with eigenvalue —kl 
solve for and /i explicitly: 



^I = 2/{Kki), 

for which the maximum criterion (6.11) reads 



s/ {nk 



') 



(6.9) 



(6.10) 



(6.11) 



we can 



(6.12) 



(6.13) 



Thus, only eigenfunction sources with kl = and kl = 2K^ are guaranteed to 
have £o > 1- The source (6.5b) has fc^ = 1^ + 2^ = 5 (with = 1), so it can 
lead to £o < as we found numerically. The same criterion (6.13) also hold for £_i to 
be bounded below by 1. 

Shaw et al [71] also derive the rigorous lower bounds 



£ > 



(6.14) 



However, these are always less than or equal to unity, so in principle they do not rule out 
the possibility than any source could be rendered inefficient by some flow. Charles R. 
Doering comments (private communication): 
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While the lower bounds (6.14) may be less than one they're greater than zero 
uniformly in Pe. That is, "unstirring" or "herding" is not something that can 
be enhanced in an unlimited manner by stirring (really a somewhat curious 
situation in my opinion!). 

The bounds (6.14) are also very permissive: for the source (6.5b) they read So ^ 1/5 
and £_i > 1/25, whereas the optimised unmixing flow in figure 6(a) has £o — -945, 
£_i ^ .642 (though the unmixing flow could be a local minimum, or could decrease the 
efliciency further at some higher Peclet number). 



7. Upper bounds on mixing efficiencies 

It is a simple matter to obtain estimates on the various mixing efficiencies. Thiffeault & 
Doering [82] used an idea of Doering & Foias [27] to obtain a simple bound on the mixing 
efficiency Zq. multiply (2.1) by an arbitrary smooth, spatially periodic 'comparison 
function' <^(cc), integrate, then integrate by parts to find 

((^(u-V + «A)(^)) = -W. (7.1) 

Here we introduced the double-bracket notation 



{F)) -.^ {F{x,t)) (7.2) 
for a space and time average, the latter defined by 

- 1 /■* 

F{x) := lim - / F{x,t')dt'. (7.3) 

t^-oo t Jo 

(We always assume that such time-averages exist.) The time derivative term from the 
advection-diffusion equation (2.1) has vanished from (7.1), since ^pdtO — dt{(p9) — 
(p{x)limt^ao9{x,t)/t — 0, since 6 is bounded. Then apply the Cauchy-Schwarz 
inequality to (7.1), to obtain 

««l>,na.WV(((-V. + .A„l (7.4) 

where the maximisation is over smooth functions (p{x). At the cost of some sharpness 
we can take the square root and then use the Minkowski inequality in the denominator, 

ieY'>max \iips}\/[i{u-Vip)Y' + K\\Aip\\). (7.5) 

We then apply Holder's inequality and find 

> max |(((^.))|/([/||V(^||^. +«:||A(^||) (7.6) 



U=(\\u\\Y (7.7) 



where 
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is proportional to the time-averaged total kinetic energy. The two bounds (7.4) and (7.6) 
both have their uses: the former is tailored to a specific velocity field, but the latter is 
a global bound valid for any stirring velocity field with bounded energy. 

To illustrate the usefulness of these estimates, we shall use (7.4) to bound the mixing 
efficiency £o, now defined to included a space-time average: where the efficiencies are 
now defined with a time average, 

£g 11^11^9 / ll^llij? ■ (7-8) 

Of course, this reduces to the earlier definition (6.1) for time-independent functions. We 
shall prove the following surprising fact mentioned in the introduction: 

An optimal way to stir a steady one- dimensional source s{x) in a periodic box, 
given a fixed time- averaged energy, is to use a spatially-uniform constant fiow 
in the x direction. 

See figure 2 for the type of source-sink distribution and fiow that we have in mind. Here, 
by optimal we mean a flow that maximises (7.8) for g = 0, or equivalently minimises 
the time- averaged variance norm {O"^)). 

Now for the proof. First take u{x, t) to be an arbitrary divergence-free vector field. 
Expand the denominator on the right in (7.4): 

{{u ■ + K^^f)) = {{u ■ V^f)) + K\{{^^f)) + 2Ki{u ■ V^)A^)). (7.9) 

Given that the source s{x) — s{x) is a function of x only, choose (p{x) = (p{x). Then 
the last term in (7.9) vanishes: 

{{{u . Vip)Aip} = iu^\xW\x))) = \iu i^p")'} = • V{ip")} = 0. (7.10) 

We also have (((u • V^pf} < (iV^Jp), where U is defined by (7.7). Hence, from (7.4) 
we have the bound 

{9'} > max {p^sf / {U' + {{A^r}) . (7.11) 

We can solve the variational problem (7.11) using its Euler-Lagrange equation, in an 
identical manner to [28,71], to find 

<(^'» > {s {i^'A^ - U^A}-' . (7.12) 

However, the right-hand side of (7.12) is exactly the variance of the periodic zero- mean 
solution to 

ue\x) = Ke"{x) + s{x), (7.13) 

that is, the steady advcction-diffusion equation for a constant flow. Hence, the constant 
flow is optimal, in the sense that any other flow with velocity norm cannot decrease 
the variance further. (This optimal solution might not be unique.) This is a surprising 
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fact: it means that any other process, even turbulence, cannot 'stir' the source-sink 
better. A constant flow also minimises the norm at fixed kinetic energy, since the 
bound for this norm is 

M%^>[s{^^'-U'^]-'s) (7.14) 

which is also saturated for the periodic zero-mean solution of (7.13). The third norm, 
, associated with £_i, is not optimised by a constant flow, reflecting the norm's 
preference for small scales. A simple bound such as (7.12) and (7.14) which does not 
depend on the details of u{x,t) cannot be derived in this case, except for particular 
classes of flows (see section 8). Nevertheless, a bound can be derived from [28,71] 

i\V-^9\^} > max(((^s))y((|Vtt- Vy' + 'u- VV<^ + /«AV(^|2», 

which after using the Minkowski inequality gives for the denominator gives what is likely 
a terrible bound: 

((I V-^^r» > max iipsf/{r' II V^ll + U || A^|| + k || VA^||)' , (7.15) 

where 7"^ = ll^'^IP is a time-averaged version of (5.7). This could in principle be 
maximised over (/?, but this is much harder than for the other norms. The important 
fact about (7.15) is that it depends on the gradient norm || Vu|| of the flow, and clearly 
the bound can be made arbitrarily small by increasing this norm. Thus flows that 
minimise this norm are likely mixing. 

If one further constrains the problem other solutions are possible (see for example 
the discussion of Plasting & Young [64] below). Nevertheless, it is surprising that the 
optimal answer in this case could be so different from a 'mixing' flow, that is, one that 
amphfles gradients of concentration (in the sense of ergodic theory - see section 4). It 
is an open question whether there exists source-sink conflgurations for which the flow 
that maximises £0 is also mixing. 

To check how sharp the bound (7.4) is for a model system, Thiffeault et al [82] 
considered the two-dimensional 'random sine flow' of Pierrehumbert [2,62]. This flow 
consists of alternating horizontal and vertical sine shear flows, with phase angles Ci(^) 
and C,2{t) G [0, 2ti] randomly chosen at each time period, r. In the flrst half of the period, 
the velocity held is 

ttW(a;,i) = V2[/(0 , sin(27rArx/L + Ci(i))) ; (7.16a) 
and in the second half-period it is 

u^'^\x,t) = V2U {sm{27rNy/ L + C2{t)) , 0), (7.16b) 

where is an integer indicating the scale of the flow. The source function used 
is s{x) = \/2sm(2-Kx/L), and wc set the integer A" = 1 for now. For the simple 
choice (fi = s, we have the efficiency bound 
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Pe 

Figure 7. Mixing efRciency £o for the random sine flow. The solid line is the upper 
bound (7.f7). The dashed line is the result of direct numerical simulations with U 
and r fixed. The dashed-dot curve plots simulation data with n and r held constant 
while varying U (after Thiffcault et al [82]). 



This bound is plotted in figure 7 against numerical simulations for the random sine 
flow. The Peclet number is varied in two ways: by varying the diffusivity n and the 
amplitude U . The bound captures the trend of the numerical simulation, especially 
as K is varied. The oscillations as U gets larger are due to the spatial periodicity of the 
domain. 

Plasting & Young [64] enhanced the bound by including the scalar dissipation rate 
as a constraint. They define the entropy production (or half the variance dissipation 
rate) as 

X:=4|V^^|^» (7.18) 

which is of course proportional to the time-average of the iJ^-norm of 9. The entropy 
production satisfies the power integral 

X=ies)). (7.19) 

Plasting & Young minimise the variance subject to both (7.1) and (7.19), taking x ^ 
given. They find the bound 

//^2vx > {mx' + HsA))is^))x+is^fis')) 

Their bound takes into account the creation of scalar gradients through the 
constraint (7.19). For the sine flow (7.16), their bound is plotted in flgure 8 in the ((6'^))- 
X plane (parabolic solid curve). The horizontal curve at the bottom is the bound (7.17). 
Notice that, if we know x, the lower bound of Plasting & Young is a vast improvement 
over(7.17). The problem is that we usually don't know x- 
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Figure 8. For the sine flow (7.16): tlie lower bound (7.20) on variance ((^?^)) (solid 
parabolic line) as a function of the dissipation The horizontal solid line is the 
bound (7.17), which assumes no knowledge of x- The vertical dashed line is the N ^ oo 
limit for x, and the dots are numerical simulation results (after Plasting & Young [64]). 



However, for the sine flow we can find x the hmit — t- oo. First, for the sine 
flow (7.16) we can compute the effective diffusivity D explicitly [34,53]: 

D = \U'^T (7.21) 

where we neglect the small molecular diffusivity. For large in (7.16), we can solve for 
a 'mean- field' temperature field [64], and obtain 

L2 2L2 ^ ^ 

for the same source s{x) = \/2sm(27fx/ L) used previously. This is the 'homogenisation 
limit,' where the scale of the source is much larger than the scale of the flow (see 
section 9). The dashed line in figure 8 shows the large N form (7.22), and the dots are 
numerical simulations by Plasting & Young for various (A^ = 1 is off scale). The 
numerical results approach the vertical line for remarkably small A^. They conjecture 
that for the sine flow as the scale separation (A^ is made larger) is increased the lower 
bound is approached. It is an open question whether the bound (7.20) can be realised 
for more general classes of flow than the single- wavenumber sine flow. 

Recently, Alexakis & Tzella [1] addressed the issue of getting bounds that reflect 
mixing instead of transport in a different way. They focused on a dissipation length 
scale id and its inverse defined by 



kl = if ■■= \\efH^/\\ef = = x/H(^l) . (7.23) 

(This length scale was denoted A in Thiffeault et al [82].) This length scale characterises 
the scale or variation of the passive scalar. They find that scale is not always 
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equivalent to the Batchelor length scale 

4 := kQU (7.24) 

where is the typical length scale of the velocity field. They introduce the ratio 

5 := 4/4 (7.25) 

where 4 is a typical length scale of the source. For example, for the sine flow (7.16) we 
have 4 ~ L/N^ and for the source s{x) = \/2 sin(27ra;/L) we have ig = L, so 5 = N~^. 
The homogenisation limit is when 5 — )■ 0, used by Plasting & Young to obtain (7.22). 
This represents the ideal of scale separation between a small-scale stirring velocity field 
and a large-scale source. (See section 9.) 

Alexakis & Tzella define the correlation ^^i^^ between the source and the 
concentration field by 

from which 



Given k and ((s^)), there are then two ways to reduce the variance: decrease ^g^s or 
decrease 4 (equivalently, increase /cd)- Decreasing is best achieved by transport, that 
is, by having a flow that rapidly carries source onto sink and vice-versa. Decreasing 4 
relies on mixing, that is, by creating small scales of the concentration fleld. Both achieve 
the same thing in the end, but in very different ways. Thus, one can target whichever 
variance- minimising method one prefers by focusing on ^ or 4- So far, our emphasis 
for the advection-diffusion problem with sources and sinks has been on decreasing 
One advantage of aiming instead to decrease 4 is that the flows obtained can be good at 
reducing the variance regardless of the precise structure of the source-sink configuration. 

A simple bound on kd was given in Thiffeault et al [82], and after being adapted 
to the two scales and Ig it reads 

klii < 5 (ci + 025 Pe-^) , Pe, := UQi^, (7.28) 

where we used a new version of the Pcclet number based on the source scale, and the 
dimensionless constants Ci and C2 depend on the shape of the source and velocity fleld. 
By working directly from the time-evolution equation for V^, Alexakis & Tzella [1] 
improve this to 

kl^l < + |^ci + 453Pe-ic2(ci + C25Pe;^). (7.29) 

where C3 is a dimensionless constant that depends on the shape of the velocity field. At 
large 5 and large Pe, this is a vast improvement over (7.28), as can be seen in figure 9. 
The constants were chosen for the sine flow with a sinusoidal source: Ci = 2\/2, C2 = 2, 
C3 = \/2 [1]. At smaller 6 the crude bound (7.28) does better, which suggests a further 
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6 

Figure 9. The upper bound (7.28) (dashed hne) and the improved bound (7.29) 
(solid hne), for various values of the Peclet number Pe„. The constants were chosen 
for the sine flow with a sinusoidal source: ci = 2-\/2, C2 = 2, C3 = ^/2. Note that 
for Pck > 1 the bounds (7.28) 'bunch up' and become independent of Pe„ (after 
Alexakis & TzeUa [1]). 

improved bound could be derived which captures both. At larger 5 the bounds (7.28) 
'bunch up' and do not improve further, whereas the bound (7.29) continues to decrease, 
achieving an asymptotic value A;j£| < C3 as Pe„ — )■ 00. 

From the upper bounds (7.28) and (7.29), Alexakis & Tzella [1] investigate the 
behaviour of as 5 varies. They identify five regimes, which are summarised in table 1. 
In regime I diffusion is very fast, so scalar gradients are set by the source scale and the 
flow is irrelevant. Regime II is a transitory regime where U and k now appear explicitly, 
but the length scale iu is still absent. Regime III is the classical Batchelor regime [12], 
where the gradients of scale as Regime IV is also Batchelor-like, in that k\ is 
proportional to f///?, but where the length scale lu has been replaced by Is-, since the 
source now has larger scales than the velocity field. It is indeed remarkable that all 
these regimes can be captured by (7.28)-(7.29). 

There is a fifth regime not captured by these bounds: this is the homogenisation 
regime when the source has much larger scale than the velocity field {5 <^ 1). We 
discussed this regime earlier in connection with the Plasting & Young bound; see also 
Majda & Kramer [53] for an extensive review, or Kramer & Keating [46] and Keating 
et al [44] for a treatment explicitly involving sources and sinks. In homogenisation 
theory the resulting effective diffusivity D usually scales as D ~ Pe"/t, where a = 2 
for shear flows (Taylor-Aris dispersion [5,78]), a = 1 for perfect chaotic mixing {D is 
then independent of k, as for the sine flow in (7.21)), and a = 1/2 for cellular flows. 
Note that, as pointed out by Alexakis & Tzella, the range of validity in 6 of regime V 
in table 1 may strictly speaking be beyond homogenisation theory: Lin et al [49] have 
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Table 1. The different regimes deduced from the bounds (7.28)-(7.29) (I-IV) and by 
homogenisation theory (V). 



regime fc^ estimate range of vaHdity note 



I 




<5>Pe 


diffusion-dominated 


II 




Pe^/^ < ^ < Pe 




III 




0(1) < (5 < Pe^/^ 


Batchelor regime 


IV 




5<0{l) 




V 




(5< min(l,Pe^-") 


homogenisation regime 



shown that 5 <^ Pe is required for the theory to apply. 
8. Dependence of norms on source-sink structure 

In this section wc discuss the results of Doering & Thiffeault [28] and Shaw, Thiffeault, 
& Doering [71], who derive bounds on the dependence of the mixing efficiencies with 
Peclet number. The 'classical' scahng for a smooth source-sink distribution is linear 
in Pe. But if the source has comphcated small-scale structures ('roughness'), then 
the efficiencies can scale anomalously with Pe, with exponents less than unity, or even 
logarithmic corrections. The specific behaviour depends on the degree of roughness, as 
characterised by the rate of decay of the power spectrum for large wavenumbers, as well 
as the dimensionality of space. 

As usual, we consider u{x,t) it to be a specified divergence-free vector field. In 
addition, we assume the following equal-time single-point statistical properties shared 
by statistically homogeneous isotropic flows (SHIFs): 



= 0, Ui{x, ■)uj{x, •) = —Sij 



IP 
d 



duj{x, •) _ Q dujjx, ■) duj{x, ■) _ 7^ 
dxk ' dxk dxu d 



(8.1) 



where overbar represents the long-time average (assumed to exist) at each point in 
space. (See [71] for a derivation of these properties.) The r.m.s. velocity U measures 
the strength of the stirring and 7 indicates the fiow field's strain or shear content. 
The ratio X = U/'y corresponds to the Taylor microscale for homogeneous isotropic 
turbulence. The Peclet number for the flow is Pe — UL/k. Note that there are flows, 
such as the 'random sine flow' (7.16), that satisfy the SHIP conditions (8.1) but are not 
genuinely isotropic [14,28,71]. Nevertheless, we will refer to flows satisfying as SHIFs 
for expediency. 

The reason for introducing SHIFs as a class of flows, in addition to their simplicity 
and physical relevance, is that the maximisation over (p in (7.4) is particularly simple. 
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for then the denominator in that equation becomes 

{{{u ■ V<p + KA<pf} = (^|V<^r + ^ (A(^)'^ . (8.2) 

Assuming a time- independent source, the simple variational problem (7.4) then 
gives [28, 71] 

£2 < A''^) 3^ 

° ~~ {s{A^-{uyK^d)A}-'s) ■ ^ ' ' 

This bound depends on the spatial structure of the source function, but not its 
amphtude; the stirring velocity field only enters through the length scale k,/U — Pe~^L. 

We can bound the small scale and large scale efficiencies £±i from (7.1) in the same 
manner after integrations by parts and application of the Cauchy-Schwarz inequality. 
For £i, 

{{ipsf = {{{uip + KVip) ■ vef <i\u^ + KVip\^}{{\ve 

so 

\2 ////.,. , .,V7,_\2 



V^|^» > max{{ips)y/{{uip + KVipf} . (8.4) 

A potentially sharper bound involving the full two-point correlation function for the 
velocity field can be obtained by formally minimising over 9 [71], but for our purposes 
the estimate (8.4) suffices. For SHIFs the denominator in (8.4) is ((/?[— A + U^]ip) and 
optimisation over (p leads to [28, 71] 

< (^(-A-)^> ^ (8.5, 

for a time-independent source. 

We can obtain a bound on £_i from (7.1) by using ^ = V • V~^9, integrating by 
parts, and using Cauchy-Schwarz: 



{{ipsf = ((V(w -Vip + KAip)- y 

< ((|Vw- Vyj + u- VV(/? + kAV<^|^))((|V" 

so that 



lfl|2\ 



((|V-'^n> >max 



V ((|Vw ■ Vv? + it ■ VV^ + kAV^P)) 
For SHIFs the denominator is {(p[-K'^A^ + {U^/d) A^ - {j'^/d) A]ip) so that 

~ {s {-A3 + (uyK^d) A2 - (7V«2rf) A}-' s) ^ ' ' 

for a time-independent source. 
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Assuming again that the fluid domain is periodic and that the source is time- 
independent, it will be helpful to rewrite the mixing efiiciency bounds (8.3), (8.5), and 
(8.6) in Fourier space: 



where we have rescaled [0,L]'^ to [0,1]'^ so that wavevector components are integer 
multiples of 27r. Now we investigate the large Peclet number behaviour of these bounds 
for a variety of classes of sources. 



8.1. Eigenfunction sources. 

The simplest class consists of sources that depend only on a single wavenumber kg, i.e., 

that are eigenfunctions of the Laplacian A with eigenvalue —k^. The bounds (8.7) then 
simplify to 

El < yfl + F^/k^ , (8.8a) 

Eo < ^/l + pJJk^ , (8.8b) 

£_i < ^ 1 + PeVA;2(i + Pe^X^k^d . (8.8c) 

Observe that each efficiency is asymptotically proportional to Pe for large Pe, 
corresponding to the expected suppression of variance if the molecular diffusivity k 
is replaced by an eddy diffusivity proportional to UL. Moreover these upper bounds are 
sharp: they may be realised by uniform flow flelds whose direction varies appropriately 
in time to satisfy the weak statistical homogeneity and isotropy conditions used in the 
analysis [64,71]. Each estimate also exhibits a decreasing dependence on the length scale 
of the source: for large Pe the bounds for the small- and intermediate-scale efficiencies 
£i and £o arc proportional to Pc/ks- 



8.2. Square-integrable sources and sinks. 

The next simplest case is when the Fourier coefficients of the source-sink distribution 
are such that the sums in the denominators of (8.7) converge in the limit as Pe — > oo. 
For example, the Fourier coefficients of smooth sources decay exponentially for large k, 
so convergence is guaranteed. We can then use the asymptotic Pe — > oo behaviour of 
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Figure 10. Plume of ash from Eyjafjallajokull volcano, Iceland, 10 May 2010 (NASA 
MODIS image). 



the mixing efficiency bounds to find 



£i<PeW ^^' :! : , (8.9a) 









k ' 


Sfcl 


2 ' 















V(A;4 + P/A2) 



These are the same Pe scahngs as above for eigenfunction sources, but the prefactors 
now depend on different combinations of length scales in the source. For instance, the 
efficiency £i depends more strongly on the high wavenumbers in the source than the 
other two efficiencies, as would be expected. 



8.3. Rough sources. 

Rough sources are very common in nature: oil gushing from an underwater well is an 
infamous recent example of a point source, as is the volcanic ash plume in figure 10. 
We can define 'rough' sources as those for which some or all the sums in the large-Pe 
bounds (8.9) are divergent. For instance, if the source is not in Li^ then the denominator 
in (8.9a) will diverge. For those cases, the Peclet number scaling may change, resulting 
in anomalous behaviour for some or all of the efficiencies (8.7). 

The extreme case is the roughest physically-meaningful sources: measure-valued 
sources such as (5-functions with nondecaying Fourier coefficients |sfc| = 0(1) as — )■ oo. 
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Then the sums in (8.9a) for £i and the denominator of (8.9b) for £o diverge in dimension 
d = 2 or 3, rendering those scahngs invahd. In this case the Pe dependence of £i 
drops out completely and all finitc-kinetic-energy stirring velocity fields are completely 
ineffective at suppressing small scale fiuctuations. To determine the high-Pe behaviour 
of £o we approximate sums by integrals. The denominator of (8.7b) is 

For d = 2 the integral in (8.10) is 

r 1^ (8 11) 

resulting in the asymptotic bound 

£o<Pe/VlogPe, d = 2. (8.12) 

Hence in dimension two there is a logarithmic correction to £o as compared to the 
square-integrable source case. 

For d — 3 the integral in (8.10) becomes 

k^dk 1 

resulting in an anomalous scaling bound 

£o<\^, d = 3. (8.14) 

This is a dramatic modification of the classical scaling. A similar analysis shows that 
the upper bound on the large scale mixing efficiency £_i ~ Pe in (8.9c) persists even 
for these roughest sources. 

We may also analyse anomalous scalings for more general rough sources where 
the Fourier spectrum \sk\ decays as k~^ with < (3 < d/2. The roughest measure- 
valued sources have /3 — while for /3 > d/2 the source is square-integrable and thus 
effectively smooth as far as these multiscale mixing efficiencies are concerned. In order 
to examine the high-Peclet-numbcr asymptotics of the bounds on the various Eg we 
estimate integrals similar to (8.10) but with an extra factor of k~'^^ in the numerator 
arising from |sfcP- The results are summarised in table 2. In = 2 the scaling for £i 
is anomalous for any degree of roughness while £o is anomalous only for the roughest 
sources with /3 = 0. In d = 3, £i is again anomalous for any degree of roughness while 
£o scales anomalously for < /3 < 1/2. For both d — 2 and 3 the bound on the large 
scale mixing efficiency £_i is always classical (i.e., linear in Pe). Of course these scalings 
neglect any large-/c cutoff for the rough sources, as discussed in the next section. 
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Table 2. Scalings of the bound on the mixing efficiency £g as functions of the source 
roughness exponent in two and three dimensions. 



fi= 2 


<l = I 


q = {) 


<1 = -1 


/3 = 


1 


Pe/(logPe)i/2 


Pe 


< /? < 1 




Pe 


Pe 


/3 = 1 


Pe/(logPe)i/2 


Pe 


Pe 


/3>1 


Pe 


Pe 


Pe 


d = 3 


/3 = 


1 


PeV2 


Pe 


< /3 < 1/2 


1 


Pe^+i/2 


Pe 


/3 = l/2 


1 


Pe/(logPe)V2 


Pe 


1/2 < /3 < 3/2 


Pc^-l/2 


Pe 


Pc 


/3 = 3/2 


Pe/(logPe)i/2 


Pe 


Pe 


/3 > 3/2 


Pe 


Pe 


Pe 



8.4- Rough sources with a cutoff. 

In nature, it can be argued that rough sources are never truly encountered: physical 
systems tend to be smooth beyond a certain small scale (as long as we stay away from 
atomic scales, but that is a different story. . . ), or at least they are modelled that way. 
With this in mind, how are the scalings derived in the previous section realised by 
sources which are only rough when seen 'from afar' but are actually smooth upon closer 
examination? Answering this will help, for instance, in understanding how such scalings 
can be observed in data for which the roughness exponent is meaningful for a limited 
range of wavenumbers. We will focus on the roughest type of sources for which Fourier 
coefficients do not decay, but the analysis is easily extended to any type of rough source 
discussed in section 8.3. 

Point-like sources of small but finite size have Fourier coefficients Sk that are 
approximately constant in magnitude up to a cutoff wavenumbcr of order 27r/£s, beyond 
which the spectrum decays as for a smooth source. We may deduce the behaviour of the 
bound on £o for such sources by inserting an upper limit at L/f^ ^ 1 into the integral 
in (8.10). For large but intermediate Peclet numbers satisfying 1 <^ Pe <C L/lg, the 
cutoff is irrelevant so the logarithmic correction (8.12) in d — 2 and the anomalous 
scahng (8.14) m. d — ?> appear. However for Pe^ L/ig, i.e., when the modified Peclet 
number based on the smallest scale in the source U£s/k,^ 1, the smooth source results 
apply and we recover the mixing efficiency bounds linear in Pe, as in (8.9). Figure (11) 
shows this scaling transition for the d = 3 case. Even in the ultimate regime where 
the source appears smooth, the prefactor in front of the high-Pe scaling bounds are 
significantly diminished by the small scales in the source: £o ^ [^og{L/£s)]~^^^ in 
d = 2, and £o < [4/L]^/2pg d = 3. 
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Pe 

Figure 11. Upper bound for the mixing efficiency £o as a function of Peciet 
number for a small source with Is = lO^^L stirred by a three-dimensional statistically 
homogeneous and isotropic flow [computed from Eq. (8.7b)]. The intermediate Pe^^^ 
scahng for K Pe < {L/ts) is evident (after Doering & Thifi:eauh [28]). 



L/50 (bound) 

- ❖ - L/50 (simulation) 

S-function (bound) 

— ^8-function (simulation) 






Figure 12. The theoretical upper bounds and simulation results for the mixing 
enhancement of a cubic source of size is = L/50, and a (5-function source. The stirring 
velocity field is a 3D version of the random sine fiow [62] (after Okabe et al [56]). 



8.5. Summary and numerical evidence. 

An important aspect of the mixing efficiency for rough source-sink distributions 
discussed above is that an anomalous large- Pe scaling implies that molecular diffusivity is 
always important, for any SHIP. To emphasise: there is no 'residual' effective diffusivity 
due to stirring in the limit of negligible molecular diffusion. Since we derived upper 
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bounds, the actual scaling could be worse (see below for numerical results). This again 
highlights a theme of this review: the source-sink distribution takes centre stage, and 
should not be treated as a secondary aspect when compared to the stirring flow. 

For actual SHIF flows, there is evidence that the upper bound scalings are somewhat 
generous for finite-size sources. Figure 12 shows results from Okabe et al [56], who used 
a 3D generalisation of Pierrehumbert's random sine flow [62] as a SHIF. The dashed 
fines are for a source of size Ig — L/50, where L is the domain size. In that case, the 
upper bound scales is a rather poor indicator of the actual flow efficiency. However, 
as the source size is made smaller the bound improves considerably: the solid lines in 
figure 12 are for a 5-function source. Chertock et al [21] have confirmed these results 
using an accurate operator-splitting method. 

Whether the results derived in this section for statistically homogeneous and 
isotropic flows (SHlFs) generalise to a wider class of flows is an important open question. 
Another open question is whether there exist SHIFs that saturate the scalings of table 2. 

9. Homogenisation theory with sources and sinks 

Homogcnisation theory is a type of multiscale analysis that exploits a large spatial scale 
separation between the stirring velocity field and the source [29,34,53,61]. We used this 
already in describing Plasting & Young's results for the sine flow (7.16) for ^ 1 in 
section 7, and discussed it in the same section in connection with regime V of Alexakis 
& Tzella. 

Here we review the results of Kramer & Keating [46] , Keating et al [44] , and Lin 
et al [49]. This series of papers has its origin in a definition of an equivalent diffusivity 
defined by Thiffeault et al [82] and generalised in Shaw et al [71]: 

D^'^'^-.^KEg (9.1) 

where Eg is just the mixing efficiency (7.8). This is the diffusivity that would be required 
in the absence of stirring to achieve the same level of suppression of the norm if^ 
with stirring. Homogenisation theory has its own effective diffusivity, D, which we 
encountered before, based on a mean-field approach and exploiting separation of scales. 
(We note that the effective diffusivity is additive on top of the molecular diffusivity, 
whereas Dq^"^^ includes the molecular diffusivity.) For large Peclet number, the effective 
diffusivity satisfies the rigorous scaling bound [6, 7, 34] D < C kFq^, whereas Dg^"^^ 
satisfies Dg^^^ < C2K,Pe for g = 1, 0, —1 (see table 2). Since both quantities claimed to 
measure essentially the same thing, both large-Pe scafings could not be right. 

The answer, of course, as explained in great detail for a model system by Lin 
et al [49], is that the homogenisation bound applies only when the concept of an 
effective homogenisation diffusivity exists, that is, in the homogenisation limit 5 <^ Pe~^, 
where the scale separation S = lu/C-s was defined in (7.25). Lin et al point out that 
the homogenisation limit 5^0 does not commute with the large Pe limit. As a 
result, the large-Pe dependence of a mixing efficiency such as £o has two distinguished 
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regimes which cross over when 5 is of order Pe. They exhibit a specific example 
where the efficiency £o ~ for fixed 5 and Pe — )■ oo, consistent with an 

upper bound hnear in Pe. For fixed Pe and 5 — ?> 0, they recover the homogenisation 
scahng £o ~ 1 + Pe^. They then introduce a modification of Batchclor's dispersion 
theory [11], caUed Dispersion-Diffusion Theory, which successfully reconciles effective 
diffusion in terms of particle dispersion and of suppression of variance of a source-sink 
distribution. 

A general treatment of homogenisation theory with sources and sinks was given 
by Kramer & Keating [46]. The starting point is a rescaled version of the advection- 
diffusion equation (2.1): we let x' — 5x, t' — SH and immediately drop the primes to 
get 

where we have assumed that the source can vary on the small spatial scale x/5 and fast 
time scale t/S'^ as well as on the larger spatial scale x and slower time scale t, whereas the 
velocity field is confined to small spatial scales. The variable x/S is assumed periodic, 
with u having space-time mean zero. The concentration 9 is expanded in the usual 
manner 

e{x, t) = e^^\t r; x,t) + 6 ^«(^, r; x, t) + e^^\t T;x,t) + ... (9.3) 

where £ := x/S, r := t/5'^ (not to be confused with the period r of the sine fiow (7.16)). 
With the fast and slow variables separated as in (9.3) we must write dt — >■ S~'^dr + dt, 
V — )■ 5~^V(_ + Va;- The advection-diffusion operator on the left-hand side of (9.2) then 
splits into three orders in 5: 

£(0) ^dr + u-V^- kA^, at order 5°; (9.4a) 
£(1) ^uW^- 2kV^ ■ V^, at order S^; (9.4b) 
£(2) ^dt- kA^, at order 5^. (9.4c) 

We have not yet posited a magnitude for the source. The simplest case that 
will yield a nontrivial self-consistent solution is to take a weak source, s — 5^s, 
where s is order one. Then at leading order we have £(0)^(0) = 0, which means 
that 9^^\^,T;x,t) = Q^^\x,t), i.e., it depends only on the slow variables. (This is 
a consequence of the solvability condition and uniqueness results — sec Lemma 3.1 
in [46].) At the next order we have = -L^^^e^^'> = -u ■ V^Q'^^l The solvability 

condition for this equation says that u • Va;©^°^ must average to zero over the small 
scales, which it does since u averages to zero and 0*^°^ does not depend on the small 
scales. Hence, we have ^^^^ = O'^^^x, t) + x- Va;©^°\ where % satisfies the so-called cell 
problem 

L^'^X^-u{tr), (9.5) 
where x has space-time mean zero. 



Using multiscale norms to quantify mixing and transport 



36 



The next and final order has = —L^^H^^^ — Z^'^^O^^^ + s, but we only require 

the solvability condition that the space-time average of the right-hand side over the 
small spatial and temporal scales vanishes. This solvability condition leads directly to 
the homogenised diffusion equation 

9te(°) = V, • {B{x, t) V,e(°)) + S{x, t) (9.6) 

where I is the unit tensor, 5* is s averaged over small space-time scales, and D is the 
tensor 

mij{x,t) := K{l+{{Vai^aj))^,r)^ S{x,t) = p))^,, (9.7) 

where the subscripts ^, r remind us that the average is over small and fast scales ^ 
and T. If the system is isotropic, the scalar effective diffusivity D we introduced earlier 
appears on the diagonal of D. 

So far everything has proceeded as one would expect: the homogenised 
equation (9.6) is exactly the standard one with an averaged source added. How can 
things go wrong and become more interesting? The most obvious way is if the small- 
scale average S in (9.7) vanishes identically. In that case we do not get a self-consistent 
equation involving the source, and we must rescale the source differently. Following 
Kramer & Keating, we set s = 6s, which is a stronger source. At order 5° nothing 
changes from before, but at order 6^ we get L^^^O'^^^ — —u • VajG^") -|- s, the source 
now making an appearance. The solvability condition is still satisfied, since the source 
averages to zero at the small scales, by assumption. We may then express the solution 
as 

9('\C r; X, t) = e«(x, t)+x- V^eW + 9,{C r; x, t) (9.8) 

where x is again the unique periodic mean-zero solution to the cell problem (9.5), 
and 9s{$,T;x,t) is the unique periodic mean-zero solution to 

£(°)^,(tr;a^,t)=s(|,r;a^,t) (9.9) 

for every {x,t). Again, this 'source cell problem' has a solution because the source 
satisfies ((s))^ ^ = 0, by assumption. The term Og gives us one extra term in the solvability 
condition for the next order, which now appears as a source instead of the vanishing 
averaged source in (9.6): 

dS^^) = V, • {B{x, t) V,e(°)) - V, • {{u 9s)) (9.10) 

where D is defined as in (9.7). Thus in this case because the small-scale mean of 
the source vanishes the source of concentration arises from the small-scale correlations 
between u and 9s that give rise to large-scale variations. To paraphrase Kramer & 
Keating, 9s is exactly the local response of the passive scalar field to the source on the 
small scale, with the large-scale variation frozen at its local value. Then u9s is the 
advective fiux of passive scalar density generated in response to the local behaviour of 
the source. It would be of great interest to generate examples of this type, and to study 
how their efficiency scales. 
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The is a third case mentioned by Kramer & Keating [46], when Va; • ^uOg])^^ 
vanishes as well as i9s}^T-- Then we must promote the strength of the source again, so 
that it comes in at order zero. Two source cell problems must then be solved for the, 
but the resulting homogenised equation looks very similar to (9.10). 

The mixing efficiencies associated with (9.7) can easily be derived [44], since in the 
absence of flow we just set D = 

El = iSA-'S}/i\{V^ ■ {B/k) ■ V^r'S\'}, (9.11a) 

£^ = iSA-JS}/i\V-' (V. • (B/k) . V^r' Sn, (9.11b) 



where now the space-time averages are over large scales. These simplify considerably 
iiB{x,t) = DI = const.: 

to = D/K, £_i=L)/k. (9.12) 

Thus, in this homogenisation limit these two efficiencies are identical, and the definition 
of equivalent diffusivity (9.1) for q = —1, is the same as the effective diffusivity D (if 
we include the additive molecular value). 

The astute reader will have noticed that we did not list £i. Indeed, Keating et 
al [44] showed that £i expressed as for (9.11) directly in terms of — is not correct 
in the context of homogenisation theory. The cause is that the gradient now takes the 
scale-separated form V = S~^'V^ + 'Vx, so the 5"^ can promote the smaller-order term 9i 
to leading other. We thus have 

we = V^^(°) + V^e(^) + 0{S^) (9.13) 

from which 

<(|V^r»,,. = «|V.^^°f + «|V,^«0>,,, + 2«V.^(°) . V,^«»,^^ + 0(5^) (9.14) 

The last term vanishes since Va;^*-^-' doesn't depend on the fast variables, and then 
using ^(0) = e(°) and e^^'> = Q'^^\x,t) + x ■ Va,e(°) we get 



in the limit as 5 — )■ 0. Inserting the steady solution 0*^°) = (Va; • (D/k) ■ Va,) ^ we 
obtain the gradient norm efficiency 

£? = {SA-''S}/{V (Va, • (D/ac) • Va,)-' S ■ (D/ac) • V (Va, • (D//^) • V,)-' S}. (9.15) 

The difference from directly trying to generalise (9.11) for the gradient norm mixing 
efficiency is that here there is an extra (D/k) sandwiched in the denominator. It is 
easier to see how this differs from (9.11) by speciahsing to D(a;,i) — Dl — const., 

£i = {D/Kf/\ (9.16) 
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and comparing to (9.12). 

Keating et al [44] suggest defining a small-scale-averaged version of the equivalent 
diffusivities, 

D^:^\x,t) ■.= Ki\VWnji\VW\^))^^^, (9.17) 

where as before 9 is the purely-diffusive solution. The advantage of definition (9.17) over 

an homogenised D is that it can be made even if there is no formal scale separation 
between the large scales and the small scales, that is, Dq^'^\x,t) characterises scalar 
dissipation due to processes at 'subgrid scales' whether or not there is a 'gap' between 
the small and large scales. In fact for g = ±1 we can go a step further and naturally 
generalise (9.17) to tensorial quantities, 

with g = — 1 being the preferred choice to relate to the effective diffusivity, because of 
the different scaling for £i in (9.16). 

10. Optimisation for the source-sink problem 

Given sources and sinks, there is an obvious optimisation problem: for fixed energy, 
which incompressible velocity field has the highest mixing efficiency £q? The stirring 
velocity field should in principle satisfy a fluid equation such as Stokes or Navier-Stokes, 
but we can also optimise over all incompressible velocity fields to get an upper bound 
on efficiency. This optimisation problem will be discussed in section 10.2. 

The presence of sources and sinks implies a different optimisation problem. If the 
position of the sources and sinks is part of the design process (as it often is for industrial 
applications), then we may try to optimise the source-sink locations as well. This is a 
less familiar problem, but one that is easier to tackle because of the structure of (2.1). 
We discuss thus this problem first in section 10.1. 

10.1. Source optimisation 

By far the easier optimisation problem is one that is less intuitive: given a stirring 
velocity field, what is the source-sink distribution which is best mixed by the flow. This 
problem was examined by Thiffeault & Pavliotis [84] . For example, suppose we have a 
room whose temperature we wish to control, and that there happens to be a predominant 
airflow in that room which is relatively unaffected by the temperature distribution. Then 
we can ask where to put heaters (sources) and windows (sinks) so that the temperature 
is as uniform as possible. 

We illustrate the optimisation procedure on the time-independent advection- 
diffusion equation, 

u{x) -Ve - kA9 = s{x), V-u = 0, (10.1) 

in f2 = [0, LY with periodic boundary conditions. The velocity field u{x) is specified. 
Both u{x) and s{x) are assumed to be sufficiently smooth. As before, we assume that 
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the source and initial condition have spatial mean zero, which implies that the scalar 
concentration also has mean zero. 

Our goal is to maximise the efficiency £q defined by (6.1), 

£J = ||(-A)'^/^^||V||(-A)^/^^|^ 

where 9 solves equation (10.1) in the absence of advection, —k,A9 = s. In maximising the 
efficiency Eg, we fix the norm of the velocity field (or equivalently, the mean kinetic 
energy of the fiow), and vary the diffusivity through the Peclet number Pe = ||tt||L/«;. 
Define the linear operators 

L :— u{x) ■ V — and L :— — kA, 

from which we can write the solution to (10.1) and to the purely-diffusive problem 
as 6* = £j~^s and 9 = L~^s, respectively. We can then rewrite the efficiency (6.1) as 

||(-A)?/2/:-is||2 (sA-'s)' ^ ' 

where the self-adjoint operators Aq and Aq are 

Aq:= L{-A)-'iL\ :=£(-A)-«£* = ft:2(-A)2-^ (10.3) 

and as before (■) denotes an average over O. To maximise we compute its variation 
with respect to s and set it equal to zero, 

= JJ^^ { (^."^ - ^^'^) = 0' (10-4) 



which implies 



A-^s^ E^A-^s. (10.5) 



This is an eigenvalue problem for the operator AqA^^. Optimal sources arc given by 
ground states of the inverse of this operator, and the normalised variance is given by 
the corresponding (first) eigenvalue. The minimisation problem has a unique minimum, 
though it may be realised by more than one source-sink distribution, in particular when 
the flow has symmetry [84] . 

The operators A~^ and A~^ are self-adjoint from L^(Q) to L'^{^); furthermore, 
they are both positive operators in L'^[fl) (restricted to functions with mean zero). 
Consequently, the generalised eigenvalue problem (10.5) has real positive eigenvalues, 
and the eigenfunctions s and s' corresponding to distinct eigenvalues arc orthogonal 
with respect to the weighted inner product {s , s') := (sA'^ ^s'). For numerical 
implementation, it is preferable to solve the equivalent self-adjoint eigenvalue problem 

{A-'/'AqA-'/') r^Elr, s If ^ ^ (io.6) 
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Pe = 0.01 Pe=l Pe=15 




Figure 13. Source distribution that optimises £o for the cellular flow (10.7), for three 
values of the Peclet numbers. Note how there is no source of heat over the stagnation 
points. The background shading shows hot (red, or dark grey) and cold (blue, or 
light grey) regions, separated by tepid regions (white). The black contour lines are 
streamlines of the flow. For both small and large Pe the optimal source converges 
to an invariant eigenf unction. In all cases there is no source of temperature over the 
elliptic stagnation points, but in the small Pe case there are sources and sinks over 
some hyperbolic points. 

for the eigenvector r, which then yields the optimal source distribution s. The advantage 
of the form (10.6) is that the self-adjoint structure of the operator is explicit. 

Our goal now is to calculate the optimal source and the corresponding mixing 
efficiency for some simple velocity fields. Notice that the operator is a diagonal 
operator in Fourier space with entries k^A;^''"^, where k = \k\ is the magnitude of 
the wavevector. For q < 2, this operator acts as a low-pass filter, suppressing high 
frequencies. 

10.1.1. Dependence on Peclet number We consider the cellular fiow on the domain VL = 
[0,L]^ = [0,27r]^ with streamfunction 

i/j{x,y) = \/2 sin X sin y (10.7) 

and velocity field u = {ux,Uy) = {dyip , —dxip) , normalised to make = 1. For this 
cellular fiow there are two independent optimal source eigenfunctions with degenerate 
optimal efficiency; the degeneracy is a consequence of the discrete symmetries of 
the fiow [84]. In figure 13 we show the source that optimises £o for three values 
of Pe. In the foreground are contour lines of the streamfunction. The optimal source 
distribution appears to become independent of Pe both for small Pe and large Pe, but 
the distributions are different. The transition between the two regimes occurs when Pe 
is of order unity. Though the two asymptotic sources are very different, they respect 
some general principles: the source is arranged for effective transport of hot onto cold 
and vice versa, and regions of high speed are favoured. In particular, note that the 
centre of the rolls has a nearly zero, fiat source distribution in all cases. 

Another perhaps surprising aspect of the small Pe solution in figure 13 is that it 
has complicated structure. In this large diffusivity limit, one would expect diffusion 
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Figure 14. For the flow with streamfunction (10.7), mixing efficiency £o ~ 1 for the 
optimal source distribution as a function of the Peclet number (after Thiffeault & 
Pavhotis [84]). 

to dominate and gradients to be smoothed out. But since our mixing efficiency (6.1) 
compares the variance to the unstirred case, which aheady has very low variance, any 
amount of improvement will count. Hence, the complicated source for small Pe in 
figure 13 only gives a minute improvement to the efficiency. The small-Pe optimal 
solution is particular in that it has some hot and cold spots localised over hyperbolic 
stagnation points. This is probably due to the high speeds along the separatrices being 
favoured, even at the cost of straddling hyperbolic stagnation points a little. 

In figure 14 we show the value of the optimal efficiency £o a^s a function of Pe. For 
large Pe, the efficiency typically scales linearly with Pe: this is the 'classical' scaling 
discussed in [28,71,82,84]. For small Pe, the optimal efficiency also converges towards 
unity linearly with Pe. 

In summary, the optimal source distribution becomes independent of Pe for both 
large and small Pe, but of course for small Pe the efficiency gain is minimal (since the Li^ 
norm of the velocity is fixed). 

10.1.2. Dependence on norm chosen Our final study will be to examine the behaviour 
of the optimal efficiency £q as q is varied in (6.1). In section 10.1.1 we used g = 0; 
now we fix Pe = 100, and allow q to vary over negative and positive values. Figure 15 
shows the optimal source distributions for q = —1, 0, and 1. For large |g| (not shown), 
the optimal source distribution converges rapidly to invariant patterns. The q = —1 
case in figure 15 (negative q) shows small, localised sources and sinks. In contrast, 
the q = I case (positive q) shows large, regular localised sources and sinks. In fact, 
what is striking about the pattern is its simplicity: it is what one might take as a guess 
at an efficient source distribution, with no added frills. Thus, a high power of q might 
be useful in situations where a simple configuration is preferable due to engineering 
constraints. The reason for the simplicity is that spatial variations in the source favour 
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Figure 15. For the flow with streamfunction (10.7), optimal source distribution s{x) 
for Pe — 100 and q = —1, 0, and 1. In all cases there are no sources or sinks of 
temperature over the stagnation points. (See the caption to figure 13 for a key to the 
background shading and contours.) 




-2 2 4 6 

q 

Figure 16. For the flow with streamfunction (10.7), mixing efficiency — 1 as a 
function of the exponent q for optimal source at Pe = 100. The optimal efhciency 
is symmetric about q = 1, and for \q\ ^ 1 it grows as (2.2) l9l (after Thiffeault & 
Pavhotis [84]). 

the diffusion operator in £, and as g — t- oo these are magnified. Thus, the source 
must remain as spatially simple as possible while trying to maximise alignment with 
the velocity. As g — i- — oo, spatial variations of the source are downplayed by the norm, 
allowing more complexity. 

Figure 16 shows how the optimal mixing efficiency varies as a function of q. 
For |g| ^ 1, the efficiency scales exponentially as 2.2l'?l. Note that the curve is symmetric 
about q = 1, which leads to a minimum there: whether this is true in general has not 
been proved, but no counterexample has been found. Thiffeault & Pavliotis [84] provide 
a partial proof by explicitly finding the symmetry between the operators ^2-5 and Ag, 
but only for small Pe. 
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10.2. Velocity field optimisation 

We turn now to a more obviously relevant problem, that of optimising the stirring 
velocity field for a given source-sink configuration. As in section 10.1, we restrict to 
the time-independent problem for simplicity. We will need the functional derivative 
of ll^lljjq, which arises from the variation 

m?m = 5\\{-^f'^^''sf = 6 {{{-Ay/^L-'sf) . (10.8) 
Since the velocity field only appears in the operator L defined in (10.1), we have 

511^11^^ = 2 {{-Ay/^L-h i-Ay/^SL-h) . (10.9) 
Using the property 6£j''^ — 5C leads to 

511^11^^ = -2 {{-Ay/^L-\s {-Ay/''L-' 6L L^^s) 

= -2{{j:-'*{-Ayj:-'s) du ■ v^-'s) , 

where we integrated by parts, used the adjoint of and substituted 5jC — 5u ■ V. 
Prom the definition (10.3) of the self-adjoint operator Aq, we can then write the 
functional derivative as 

l^-^--{A-h)VL~^s. (10.10) 

To formulate an optimisation problem, we also need to add constraints on it: 
incompressibility and fixed energy. This is done in the usual manner by considering 
Lagrange multipliers in the extended functional 

= m\%. + ¥ i\M" -u') + {^y-u). (10.11) 

Here // and i'{x) are the Lagrange multipliers, with v a function of space since V-it = is 
a pointwise constraint. The functional derivative of (^[w] then gives the Euler-Lagrange 
equation, 

- - {A-\s) VL-h + nu-Wu = 0. (10.12) 
ou ' 

to be solved for u for given s. Note that (10.12) is profoundly nonlinear in u, since it 
enters in the nonlocal operators L^^ and A~^. 

Let us restrict to the two-dimensional case, where we can introduce a 
streamfunction t/j with u — {dyijj, —dxi/j). Then taking the curl of (10.12) yields 

fxAiJj+ [A-^s , L'h] = 0, (10.13) 

where 

[f , g] = dj dyg - dyf d^g. (10.14) 

Equation (10.13) is a nonlinear eigenvalue problem that can be solved in several ways. 
A direct approach is to start with an initial guess for ip and and compute the residual 
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(a) (b) 

Figure 17. (a) Optimal stirring velocity field (streamlines) for the source cos a; cosy, 
for Pe = 10. (b) Dependence on Peclet number of the optimal mixing efficiency fig- 
For small Pe the optimal streamfunction approaches \/2 sin a; sin y. The dashed line is 
the upper bound (10.15). 



vector, that is the amount by which (10.13) fails to be satisfied. We also append the 
constraint (11^11^ — U"^) to the residual vector. Then use a multidimensional nonlinear 
solver such as Matlab's f solve that finds zeroes of the residual vector, by adjusting the 
vector {ip, fi). Here has been discretised in some way, either by specifying it on a grid 
or by expanding it Fourier series. 

Figure 17(a) shows a solution of (10.13) with q = 0, U = 1, and Pe = 10 for the 
source cos a; cos |/. The efficiency corresponding to this solution is £o = 1.46. The fiow 
is close to the standard cellular fiow (10.7), but with a fiattened core where velocities 
are smaller. Whether this is truly optimal is a difficult question: (10.13) has many 
solutions with different /i, and unlike the source optimisation case (which leads to a 
linear eigenvalue problem) there is no simple way of finding minimising solutions (but 
see the upper bound (10.15) below). It is an open challenge to characterise the solutions 
of (10.12) and (10.13) more thoroughly. The solution in figure 17(a) was obtained from 
the initial guess a/2 sinx siny, and all other initial conditions examined gave larger 
values of the norm. 

Figure 17(b) shows the dependence on Pe of the optimal mixing efficiency £o, for 
the fixed source-sink distribution sinxsiny of figure 17(a). For smaller values of Pe, 
the optimal £o ~ 1 is proportional to Pe^: this is the diffusion-dominated regime, where 
stirring only has a small effect. The optimal solution converges to a/2 sinxsiny for 
small Pe. For larger values of Pe, the optimal solution recovers the 'classical' upper 
bound scaling, linear in Pe. 

It is instructive to compare the optimal efficiency plotted in figure 17(b) with 
the 'global bound' (7.6). This requires a choice of comparison function and the 
simplest is to take f{x) = s{x) = cos x cosy. We then have {ips) = 1/4, ||Vv5||^oo = 1 




Figure 18. Left: a sphere heated at the equator and cooled at the poles. Right: 
longitude-latitude plot of the streamlines of a steady flow that maximises £o, along 
with the temperature field in the background (from O'Rourke [58]). 



and ||A(y9|| = 1 in (7.6). After normalising by the purely-diffusive solution we obtain 

U Pe 

£o-l<- = 7^, (hereL = 27r), (10.15) 

which is valid even if the velocity field is allowed to be time-dependent, in which case U 
is defined as in (7.7). The bound (10.15) is plotted as a dashed line in figure 17(b), 
where we see that our optimal solution is remarkably close to the global upper bound 
for large Pe (about 37.5% above the optimal solution). This shows that the series 
of inequalities required to obtain (10.15) do not cause too much loss of sharpness for 
large Pe, but it also implies that even allowing for arbitrary time- dependent stirring 
cannot improve £o by very much. The bound (7.6) thus helps to determine if our local 
optimal solution is anywhere close to being a global optimum. 



11. Discussion 



The literature on stirring, mixing, and transport is enormous, and of course we only 
covered a small corner of it, focusing on direct uses of norms. The decay of variance 
itself is an object of study, and predicting its decay rate in terms of flow characteristics 
has been a long-term goal of the theory of mixing. The dominant approaches are the 
local theory, based on dynamical systems quantities such as the statistics of finite-time 
Lyapunov exponents [2,3,8,33,72,79,88], and the global theory, which requires a more 
thorough analysis of the advection-diffusion equation [35,38,41,52,62,63,65,66,77,80,93]. 
Neither of these approaches is particularly well-suited to optimisation, so we do not 
discuss them here. It is also important to note that there are many other measures 
of mixing beyond norms — see for instance [13,22,37,39,47,60,74,75,83,89,90] and 
references therein. 

We have reviewed the reasons why various norms were desirable for studying mixing, 
and how they could be used to optimise the rate of decay for the initial-value, freely- 
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decaying problem. We then focused for most of the review on the use of norms in 
the presence of sources and sinks in the long-time limit where the system attain an 
equilibrium. Our examples used simple periodic geometries. 

Let us briefly discuss the work of Constantin et al [22], which contains some of 
the most important recent rigorous results (we leave out a few technical conditions). 
Their focus is on flows that are relaxation- enhancing. These are steady flows with the 
property that for any given t > and £ > 0, it is possible to increase the amplitude 
of stirring to make < e. (As usual we assume mean-zero functions.) In other 

words, it is at least possible to achieve an arbitrary level of mixing, measured according 
to the norm, by stirring hard enough. Weakly mixing incompressible flows are 
always relaxation-enhancing. Constantin et al prove that an incompressible velocity 
field u{x) is relaxation-enhancing if and only if tt ■ V has no eigenfunctions in 
other than the constant function. Indeed, if there exists such eigenfunctions then most 
initial conditions will contain some admixture of them, which will not decay. Current 
examples of relaxation-enhancing flows are, however, not very physical. It is possible 
that extending the work of Constantin et al to time-dependent velocity flelds would 
greatly increase the range of flows that are relaxation-enhancing, but this is likely to be 
difficult. 

There are a number of open problems and areas for further study. We mention a 

few: 

• Find the optimal velocities for more complicated flows, and reflne the numerical 
methods needed to do so; for example O'Rourke has recently examined optimal 
transport flows on the sphere [58] (see figure 18). 

• Understand better the transition between 'transporting' and 'mixing' fiows. The 
fiows that minimise the norms in the presence of sources and sinks are very different 
from mixing fiows. Is this a fiaw in the measure? What would be a better measure? 
Are there fiows and source-sink distributions which are both optimal in the sense 
of minimising norms, but are also optimally mixing in the sense of ergodic theory? 

• Since all norms for g < act as 'mix-norms,' that is, they decay if the 
flow is mixing, which one is best? Do we need to select q according to particular 
applications? Mathew et al [54,55] used q — —1/2, and in this review we focused 
more on q — —1, but there is no clear reason to choose either one at this point. On 
the plus side, however, both work quite well. 

• In a similar vein as the previous problem, it is not known whether there are 

any advantages in using even more general Sobolev norms on the space W^''^, 
1 < p < oo, rather than on W^'"^ = H'^. Values of p other than 2 are rarely used in 
the context of mixing, except for p = oo which is common (see for example [90]). 
Some rigorous results, such as [22] , do not depend on p, which suggests the choice 
of p matters little. 

• Do the results on 'roughness' of the source-sink distribution carry over to open 
flows, where the sources and sinks can be regarded as distributed on the boundary 
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rather than in the bulk [81]? 

• There has been some work on using the norm approach to quantify mixing in more 
complex systems with reactions, for example a decaying passive scalar [70], the 
Fisher equation [14], and the Cahn-Hilliard equation [57]. In this case how do 
reaction rates, etc., depend on the source-sink structure? What kinds of flows 
optimise reaction rates? 
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Appendix A. Proof of mix-norm theorem 

In this appendix we prove the theorem presented at the end of section 4 relating weak 
convergence to the negative Sobolev norms. The proof is from Lin et al [50]. Write the 
norm for H'^{Q) as 



where = (1 + k^L^ for the norm (3.1). Suppose that f{-,t) is uniformly bounded 
in L^{^), so that ||/(-,OII ^ limt^oo !!/(•) ^) ^ for some ^ < 0. Then for 

any g e L'^(n), 




(A.l) 



k 




Given e > 0, first choose K{e) such that 




(A.2) 
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then choose ^(e) such that 

\k<K{e) ) 

We then have 

|(/,5)|<i(l + C-M|/||)e<e, t>T{e\ 

which imphes that / converges weakly to zero as i — > oo. (This is true even for = 0.) 

Conversely, suppose ||/(-, i)|| < C for all t and limt_^oo {f : g) ^ ^ for all g G L^(Q). 
By choosing g — exp(— ifc-cc) we see that all the Fourier coefficients fk{t) — >■ as i — >■ oo. 
Also, because \\f{;t)\f = Efel/feWP < then each \fk{t)\ < C for all t. 

We have 

<E^?^i^i'+^^^ii/ii'- (A-4) 

For any e > 0, we can choose K{e) such that X^^!^ \\f\\ < X'j^^C < e/2 for k > K{e) (this 
requires q < 0). For any finite K, Efc<^ -^i'^'*l/fe(^)P — > as i — > oo, so there exists ^(e) 
such that Efc<if(e) ^?l/fe(^)P < e/2, for all t > T{e). From (A.4) we obtain \\f\\%, < e 
for alH > ^(e), which proves the result. 
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